# -*- coding: utf-8 -*-
"""
涡旋轨迹追踪核心程序
功能描述：基于NC格式的涡旋检测数据（含位置、振幅、有效半径、轮廓等信息），实现涡旋轨迹的连续追踪，
         支持分裂/合并事件检测、轨迹中断续接（宽容期机制），并提供轨迹统计、可视化与结果保存功能。
适用数据：NC格式涡旋文件（需包含time/latitude/longitude/amplitude/effective_radius等基础变量，
         可选effective_contour_latitude/longitude轮廓变量）
依赖库：os/re（文件路径处理）、numpy/pandas（数值计算）、h5netcdf（NC文件读写）、matplotlib/Basemap（可视化）、
        logging（日志）、datetime（时间处理）、collections（数据结构）、tqdm（进度条）、numba（加速计算）、gdal（地理数据）
"""

# --------------------------
# 1. 导入依赖库
# --------------------------
import os
import re
import numpy as np
import pandas as pd
import h5netcdf
import matplotlib.pyplot as plt
from matplotlib.lines import Line2D
import matplotlib.patches as patches
import logging
from datetime import datetime
from collections import defaultdict, Counter
from mpl_toolkits.basemap import Basemap
from tqdm import tqdm
from numba import njit, prange  # 用于加速数值计算（如相似度矩阵）
from osgeo import gdal, osr  # 预留地理坐标转换接口（当前用于经纬度校验）

# 配置matplotlib中文字体（解决中文显示乱码问题）
plt.rcParams["font.family"] = ["SimHei", "WenQuanYi Micro Hei", "Heiti TC", "Arial Unicode MS"]
plt.rcParams['axes.unicode_minus'] = False  # 正确显示负号（避免负号显示为方块）

# --------------------------
# 2. 全局参数配置（可根据数据场景调整）
# --------------------------
max_daily_distance = 100  # 涡旋日均最大移动距离约束（单位：km），增大以适配长距离移动涡旋（原500km→100km为示例调整，需按实际数据优化）
base_threshold = 0.48     # 涡旋直接匹配阈值（相似度≥该值视为同一涡旋），降低阈值可提高匹配灵敏度（原0.48为基础值，需按数据精度调整）
split_threshold = 0.3     # 涡旋分裂事件阈值（单个前序涡旋匹配多个当前涡旋时的相似度下限）
merge_threshold = 0.3     # 涡旋合并事件阈值（多个前序涡旋匹配单个当前涡旋时的相似度下限）
tolerance_steps = 1       # 轨迹中断宽容期（允许中断的最大时间步数，超过则视为轨迹终止，用于处理临时检测缺失）

# --------------------------
# 3. 日志系统配置（用于追踪程序运行状态与错误）
# --------------------------
def setup_logger():
    """
    初始化日志系统，同时输出至控制台与文件（eddy_tracking.log）
    日志级别：INFO（常规运行信息）、WARNING（潜在问题）、ERROR（错误信息）、DEBUG（调试细节）
    返回：配置完成的logger对象
    """
    logger = logging.getLogger('eddy_tracker')  # 日志器名称（唯一标识）
    logger.setLevel(logging.INFO)  # 全局日志级别（INFO及以上级别会被记录）
    logger.propagate = False  # 禁止日志重复传播（避免与其他库日志冲突）

    # 定义日志格式（包含时间、级别、信息）
    formatter = logging.Formatter('%(asctime)s - %(levelname)s - %(message)s')

    # 1. 文件处理器（将日志写入文件）
    file_handler = logging.FileHandler('eddy_tracking.log', encoding='utf-8')
    file_handler.setFormatter(formatter)
    file_handler.setLevel(logging.DEBUG)  # 文件日志记录DEBUG及以上（包含更详细调试信息）

    # 2. 控制台处理器（将日志输出至控制台）
    console_handler = logging.StreamHandler()
    console_handler.setFormatter(formatter)
    console_handler.setLevel(logging.INFO)  # 控制台仅显示INFO及以上（避免调试信息冗余）

    # 绑定处理器到日志器
    logger.addHandler(file_handler)
    logger.addHandler(console_handler)

    return logger

# 初始化日志器（全局唯一）
logger = setup_logger()

# --------------------------
# 4. 数据读取与预处理模块（加载NC涡旋数据并标准化格式）
# --------------------------
def read_eddy_file(file_path):
    """
    读取单个NC格式涡旋文件，提取涡旋核心参数与轮廓信息（若存在）
    参数：
        file_path (str): NC文件路径（如"D:/EddyData/2023-01-01_eddy.nc"）
    返回：
        list[dict] | None: 涡旋列表（每个元素为涡旋字典，含基础参数与轮廓信息）；读取失败返回None
    异常：
        捕获文件读取、变量缺失等异常，记录ERROR日志
    """
    try:
        # 以只读模式打开NC文件（h5netcdf支持HDF5格式NC文件）
        with h5netcdf.File(file_path, 'r') as ds:
            # 获取涡旋数量（NC文件中"obs"维度通常表示观测样本数，即涡旋个数）
            obs_count = len(ds.dimensions.get('obs', []))
            eddies = []  # 存储当前文件的所有涡旋

            # 读取基础涡旋参数（必选变量，若缺失会抛出KeyError，需确保NC文件格式规范）
            times = ds['time'][:]  # 时间戳（通常为Unix时间或儒略日）
            latitudes = ds['latitude'][:]  # 涡旋中心点纬度（单位：度）
            longitudes = ds['longitude'][:]  # 涡旋中心点经度（单位：度）
            amplitudes = ds['amplitude'][:]  # 涡旋振幅（单位：m，表征涡旋强度）
            effective_radii = ds['effective_radius'][:]  # 有效半径（单位：km，表征涡旋大小）
            effective_areas = ds['effective_area'][:]  # 有效面积（单位：km²，可选变量，部分文件可能命名不同）

            # 读取轮廓数据（可选变量，用于更精确的涡旋形态匹配）
            has_contour = 'effective_contour_latitude' in ds.variables
            if has_contour:
                # 轮廓经纬度通常存储为整数（如放大100倍以保留精度），需转换为度
                contour_lats = ds['effective_contour_latitude'][:] / 100
                contour_lons = ds['effective_contour_longitude'][:] / 100
                # 计算轮廓经纬度均值（替代中心点，提升匹配稳定性，避免单点误差）
                contour_lat_means = np.mean(contour_lats, axis=1)
                contour_lon_means = np.mean(contour_lons, axis=1)
            else:
                # 无轮廓数据时，使用中心点经纬度作为替代
                contour_lats = None
                contour_lons = None
                contour_lat_means = latitudes
                contour_lon_means = longitudes

            # 从文件名提取日期（假设文件名含"YYYY-MM-DD"格式日期，如"eddy_2023-01-01.nc"）
            date_match = re.search(r'(\d{4}-\d{2}-\d{2})', file_path)
            date_str = date_match.group(1) if date_match else 'unknown'  # 无日期则标记为unknown

            # 逐个处理每个涡旋，构建标准化字典
            for idx in range(obs_count):
                eddy = {
                    'obs_id': idx,  # 当前文件内的涡旋索引（本地ID）
                    'time': times[idx],  # 时间戳
                    'latitude': latitudes[idx],  # 中心点纬度
                    'longitude': longitudes[idx],  # 中心点经度
                    'amplitude': amplitudes[idx],  # 振幅
                    'effective_radius': effective_radii[idx],  # 有效半径
                    'effective_area': effective_areas[idx],  # 有效面积
                    'contour_lat_mean': contour_lat_means[idx],  # 轮廓纬度均值
                    'contour_lon_mean': contour_lon_means[idx],  # 轮廓经度均值
                    'date': date_str,  # 日期字符串（YYYY-MM-DD）
                    'contour_latitude': None,  # 完整轮廓纬度（列表）
                    'contour_longitude': None  # 完整轮廓经度（列表）
                }

                # 若存在轮廓数据，补充完整轮廓信息
                if has_contour:
                    eddy['contour_latitude'] = contour_lats[idx]
                    eddy['contour_longitude'] = contour_lons[idx]

                eddies.append(eddy)

            # 记录DEBUG日志（仅调试时可见，说明文件读取结果）
            logger.debug(f"成功读取文件 {os.path.basename(file_path)}，包含 {obs_count} 个涡旋")
            return eddies

    except Exception as e:
        # 捕获所有异常（如文件不存在、变量缺失、格式错误），记录错误日志
        logger.error(f"读取文件 {file_path} 失败: {str(e)}", exc_info=True)  # exc_info=True记录堆栈信息，便于调试
        return None


def load_eddy_series(file_dir):
    """
    加载指定目录下所有NC涡旋文件，按日期排序构建时间序列涡旋数据
    参数：
        file_dir (str): 涡旋NC文件目录（如"D:/EddyData"）
    返回：
        list[list[dict]]: 时间序列涡旋数据（外层列表按时间步排序，内层列表为对应时间步的涡旋列表）
    说明：
        1. 仅处理后缀为".nc"的文件；2. 按文件名中的日期（YYYYMMDD格式）排序，确保时间连续性
    """
    # 筛选目录下所有NC文件
    eddy_files = [f for f in os.listdir(file_dir) if f.lower().endswith('.nc')]
    if not eddy_files:
        logger.warning(f"目录 {file_dir} 下未找到NC文件")
        return []

    # 定义日期提取函数（用于文件排序）：从文件名提取YYYYMMDD格式日期
    def extract_date(filename):
        match = re.search(r'(\d{4})(\d{2})(\d{2})', filename)
        if match:
            # 转换为datetime对象，便于排序
            return datetime(int(match.group(1)), int(match.group(2)), int(match.group(3)))
        return datetime.min  # 无日期文件排在最前

    # 按日期升序排列文件（确保时间步连续）
    eddy_files.sort(key=extract_date)
    eddy_series = []  # 存储时间序列涡旋数据

    # 逐个加载文件并添加时间步信息
    for time_step, filename in enumerate(eddy_files):
        file_path = os.path.join(file_dir, filename)
        eddies = read_eddy_file(file_path)
        if eddies is None:
            continue  # 跳过读取失败的文件

        # 为每个涡旋添加时间步索引（全局时间步，从0开始）
        for eddy in eddies:
            eddy['time_step'] = time_step

        eddy_series.append(eddies)
        # 记录INFO日志（常规运行信息，说明每个文件的加载结果）
        logger.info(f"加载时间步 {time_step} 文件: {filename}，包含 {len(eddies)} 个涡旋")

    return eddy_series

# --------------------------
# 5. 经纬度转换工具（处理经度范围不一致问题）
# --------------------------
def convert_lon_360_to_180(lon):
    """
    将360°范围经度（0~360°）转换为180°范围（-180~180°），适配Basemap等可视化库
    参数：
        lon (np.ndarray | float): 360°范围经度（单个值或数组）
    返回：
        np.ndarray | float: 180°范围经度
    示例：
        convert_lon_360_to_180(190) → -170；convert_lon_360_to_180(350) → -10
    """
    return ((lon + 180) % 360) - 180


def check_and_convert_lon(eddy_series):
    """
    检查时间序列涡旋数据的经度范围，若为360°范围则转换为180°范围
    参数：
        eddy_series (list[list[dict]]): 时间序列涡旋数据（load_eddy_series的输出）
    返回：
        list[list[dict]]: 经度转换后的时间序列涡旋数据（若无需转换则返回原数据）
    说明：
        1. 检查所有涡旋的中心点经度、轮廓均值经度、完整轮廓经度；
        2. 若存在<-180°或>180°的经度，触发转换（避免可视化时经度跨接问题）
    """
    # 收集所有经度数据（用于判断范围）
    all_lons = []
    for eddies in eddy_series:
        for eddy in eddies:
            all_lons.append(eddy['longitude'])
            all_lons.append(eddy['contour_lon_mean'])
            # 若存在完整轮廓，添加所有轮廓经度
            if eddy['contour_longitude'] is not None:
                all_lons.extend(eddy['contour_longitude'])

    # 无经度数据时直接返回（异常情况，如数据为空）
    if not all_lons:
        logger.warning("无经度数据可检查，直接返回原数据")
        return eddy_series

    # 计算经度范围
    min_lon = min(all_lons)
    max_lon = max(all_lons)

    # 判断是否需要转换（360°范围特征：max_lon>180°且min_lon>=0°）
    if min_lon < -180 or max_lon > 180:
        logger.info(f"检测到经度范围: [{min_lon:.2f}, {max_lon:.2f}]，转换为[-180, 180]范围")
        converted_series = []
        for eddies in eddy_series:
            converted_eddies = []
            for eddy in eddies:
                # 转换中心点经度
                eddy['longitude'] = convert_lon_360_to_180(eddy['longitude'])
                # 转换轮廓均值经度
                eddy['contour_lon_mean'] = convert_lon_360_to_180(eddy['contour_lon_mean'])
                # 转换完整轮廓经度（若存在）
                if eddy['contour_longitude'] is not None:
                    eddy['contour_longitude'] = convert_lon_360_to_180(eddy['contour_longitude'])
                converted_eddies.append(eddy)
            converted_series.append(converted_eddies)
        return converted_series

    # 无需转换时记录信息
    logger.info(f"经度范围: [{min_lon:.2f}, {max_lon:.2f}]，符合[-180, 180]范围，无需转换")
    return eddy_series

# --------------------------
# 6. 相似度计算模块（Numba加速，核心匹配逻辑）
# --------------------------
@njit(fastmath=True, cache=True)
def haversine_distance(lat1, lon1, lat2, lon2):
    """
    计算两点间的哈夫曼距离（球面上的最短距离，适配地球曲率）
    参数：
        lat1 (float): 点1纬度（度）
        lon1 (float): 点1经度（度）
        lat2 (float): 点2纬度（度）
        lon2 (float): 点2经度（度）
    返回：
        float: 两点间距离（单位：km）
    加速：
        @njit装饰器（Numba）将函数编译为机器码，比纯Python快10~100倍；
        fastmath=True启用浮点计算优化，cache=True缓存编译结果
    """
    R = 6371.0  # 地球平均半径（单位：km）
    # 转换为弧度（三角函数计算需弧度输入）
    lat1, lon1, lat2, lon2 = np.radians(lat1), np.radians(lon1), np.radians(lat2), np.radians(lon2)

    # 哈夫曼公式计算距离
    dlat = lat2 - lat1  # 纬度差
    dlon = lon2 - lon1  # 经度差
    a = np.sin(dlat / 2) ** 2 + np.cos(lat1) * np.cos(lat2) * np.sin(dlon / 2) ** 2
    c = 2 * np.arctan2(np.sqrt(a), np.sqrt(1 - a))  # 圆心角
    return R * c  # 距离 = 半径 × 圆心角


@njit(parallel=True, fastmath=True, cache=True)
def calculate_similarity_matrix(lats1, lons1, radii1, amps1, areas1,
                                lats2, lons2, radii2, amps2, areas2,
                                max_daily_distance=50):
    """
    计算两个时间步涡旋间的相似度矩阵（核心匹配逻辑）
    参数：
        lats1/lons1 (np.ndarray): 前序时间步涡旋的纬度/经度（轮廓均值，长度=N）
        radii1/amps1/areas1 (np.ndarray): 前序时间步涡旋的半径/振幅/面积（长度=N）
        lats2/lons2 (np.ndarray): 当前时间步涡旋的纬度/经度（轮廓均值，长度=M）
        radii2/amps2/areas2 (np.ndarray): 当前时间步涡旋的半径/振幅/面积（长度=M）
        max_daily_distance (float): 涡旋日均最大移动距离（超过则相似度为0，默认50km）
    返回：
        np.ndarray: 相似度矩阵（形状[N, M]，值范围0~1，越大越相似）
    加速：
        parallel=True启用多线程并行计算（prange替代range实现并行循环），适合大规模数据
    相似度计算逻辑：
        总相似度 = 2.0×距离得分 + 0.3×振幅相似度 + 0.3×面积相似度（权重可按涡旋特性调整）
        1. 距离得分：基于哈夫曼距离与涡旋半径的比值，用指数衰减（减少距离对相似度的过度惩罚）；
        2. 振幅/面积相似度：取两者最小值与最大值的比值（表征大小/强度一致性）
    """
    n_prev = len(lats1)  # 前序时间步涡旋数
    n_curr = len(lats2)  # 当前时间步涡旋数
    similarity = np.zeros((n_prev, n_curr))  # 初始化相似度矩阵

    # 并行循环计算每个前序涡旋与当前涡旋的相似度
    for i in prange(n_prev):
        for j in range(n_curr):
            # 1. 计算哈夫曼距离，超过最大移动距离则相似度为0
            dist = haversine_distance(lats1[i], lons1[i], lats2[j], lons2[j])
            if dist > max_daily_distance:
                similarity[i, j] = 0.0
                continue

            # 2. 计算距离得分（放宽距离衰减系数至2.5，原2.2，提升长距离匹配灵敏度）
            avg_radius = (radii1[i] + radii2[j]) / 2  # 平均半径（避免单个半径为0的情况）
            if avg_radius > 0:
                dist_score = np.exp(-(dist / (2.5 * avg_radius)) ** 2)
            else:
                dist_score = 0.0  # 半径为0时距离得分无效

            # 3. 计算振幅相似度（避免除以0，振幅为0时相似度为0）
            amp_max = max(amps1[i], amps2[j])
            amp_ratio = min(amps1[i], amps2[j]) / amp_max if amp_max > 0 else 0.0

            # 4. 计算面积相似度（同理，避免除以0）
            area_max = max(areas1[i], areas2[i])
            area_ratio = min(areas1[i], areas2[j]) / area_max if area_max > 0 else 0.0

            # 5. 计算总相似度（权重：距离占比最高，其次是振幅和面积，总和>1便于后续阈值调整）
            similarity[i, j] = 2.0 * dist_score + 0.3 * amp_ratio + 0.3 * area_ratio

    return similarity

# --------------------------
# 7. 涡旋匹配与轨迹追踪模块（核心业务逻辑）
# --------------------------
def match_eddies_with_split_merge(prev_eddies, curr_eddies, threshold=base_threshold,
                                  max_daily_distance=max_daily_distance,
                                  split_threshold=split_threshold, merge_threshold=merge_threshold):
    """
    单时间步涡旋匹配（支持直接匹配、分裂、合并事件检测）
    参数：
        prev_eddies (list[dict]): 前序时间步涡旋列表
        curr_eddies (list[dict]): 当前时间步涡旋列表
        threshold (float): 直接匹配阈值（默认使用全局base_threshold）
        max_daily_distance (float): 日均最大移动距离（默认使用全局max_daily_distance）
        split_threshold (float): 分裂事件阈值（默认使用全局split_threshold）
        merge_threshold (float): 合并事件阈值（默认使用全局merge_threshold）
    返回：
        list[tuple]: 匹配结果列表（每个元素为(i,j,type)，i=前序涡旋索引，j=当前涡旋索引，type=匹配类型）
        匹配类型：'direct'（直接匹配）、'split'（分裂）、'merge'（合并）
    匹配逻辑：
        1. 先按相似度降序进行一对一直接匹配（避免多对多冲突）；
        2. 未匹配前序涡旋若匹配多个当前涡旋（≥2个），判定为分裂；
        3. 未匹配当前涡旋若匹配多个前序涡旋（≥2个），判定为合并；
    """
    # 若任一时间步无涡旋，直接返回空匹配结果
    if len(prev_eddies) == 0 or len(curr_eddies) == 0:
        logger.debug("前序或当前时间步无涡旋，无匹配结果")
        return []

    n_prev, n_curr = len(prev_eddies), len(curr_eddies)

    # 提取前序时间步涡旋特征（用于相似度计算）
    prev_lats = np.array([e['contour_lat_mean'] for e in prev_eddies])
    prev_lons = np.array([e['contour_lon_mean'] for e in prev_eddies])
    prev_radii = np.array([e['effective_radius'] for e in prev_eddies])
    prev_amps = np.array([e['amplitude'] for e in prev_eddies])
    prev_areas = np.array([e['effective_area'] for e in prev_eddies])

    # 提取当前时间步涡旋特征
    curr_lats = np.array([e['contour_lat_mean'] for e in curr_eddies])
    curr_lons = np.array([e['contour_lon_mean'] for e in curr_eddies])
    curr_radii = np.array([e['effective_radius'] for e in curr_eddies])
    curr_amps = np.array([e['amplitude'] for e in curr_eddies])
    curr_areas = np.array([e['effective_area'] for e in curr_eddies])

    # 计算相似度矩阵（调用Numba加速函数）
    similarity_matrix = calculate_similarity_matrix(
        prev_lats, prev_lons, prev_radii, prev_amps, prev_areas,
        curr_lats, curr_lons, curr_radii, curr_amps, curr_areas,
        max_daily_distance=max_daily_distance
    )

    # 初始化匹配结果与已使用标记（避免重复匹配）
    matches = []
    used_prev = set()  # 已匹配的前序涡旋索引
    used_curr = set()  # 已匹配的当前涡旋索引

    # 将相似度矩阵展平并按降序排序（获取全局相似度最高的匹配对）
    flat_indices = np.argsort(similarity_matrix.flatten())[::-1]

    # 1. 一对一直接匹配（优先匹配相似度最高的非冲突对）
    for idx in flat_indices:
        # 还原矩阵索引（展平索引→二维索引）
        i = idx // n_curr  # 前序涡旋索引
        j = idx % n_curr   # 当前涡旋索引

        # 若两者均未使用且相似度超过阈值，判定为直接匹配
        if i not in used_prev and j not in used_curr and similarity_matrix[i, j] > threshold:
            # 二次校验距离（避免相似度计算中可能的误差）
            dist = haversine_distance(prev_lats[i], prev_lons[i], curr_lats[j], curr_lons[j])
            if dist <= max_daily_distance:
                matches.append((i, j, 'direct'))
                used_prev.add(i)
                used_curr.add(j)

    # 2. 分裂事件检测（前序1个→当前多个）
    for i in range(n_prev):
        if i in used_prev:
            continue  # 跳过已匹配的前序涡旋

        # 收集当前时间步中相似度超过分裂阈值的未匹配涡旋
        candidates = []
        for j in range(n_curr):
            if j not in used_curr and similarity_matrix[i, j] > split_threshold:
                dist = haversine_distance(prev_lats[i], prev_lons[i], curr_lats[j], curr_lons[j])
                if dist <= max_daily_distance:
                    candidates.append((j, similarity_matrix[i, j]))

        # 若候选数≥2，判定为分裂（取前3个相似度最高的候选，避免过多分支）
        if len(candidates) >= 2:
            candidates.sort(key=lambda x: x[1], reverse=True)  # 按相似度降序
            for j, _ in candidates[:3]:
                matches.append((i, j, 'split'))
                used_prev.add(i)
                used_curr.add(j)
            logger.info(f"检测到分裂事件: 前序涡旋 {i} → 当前涡旋 {[j for j, _ in candidates[:3]]}")

    # 3. 合并事件检测（前序多个→当前1个）
    for j in range(n_curr):
        if j in used_curr:
            continue  # 跳过已匹配的当前涡旋

        # 收集前序时间步中相似度超过合并阈值的未匹配涡旋
        candidates = []
        for i in range(n_prev):
            if i not in used_prev and similarity_matrix[i, j] > merge_threshold:
                dist = haversine_distance(prev_lats[i], prev_lons[i], curr_lats[j], curr_lons[j])
                if dist <= max_daily_distance:
                    candidates.append((i, similarity_matrix[i, j]))

        # 若候选数≥2，判定为合并（取前3个相似度最高的候选）
        if len(candidates) >= 2:
            candidates.sort(key=lambda x: x[1], reverse=True)
            for i, _ in candidates[:3]:
                matches.append((i, j, 'merge'))
                used_prev.add(i)
                used_curr.add(j)
            logger.info(f"检测到合并事件: 前序涡旋 {[i for i, _ in candidates[:3]]} → 当前涡旋 {j}")

    # 记录匹配统计信息
    logger.info(
        f"单时间步匹配结果: 成功匹配 {len(matches)} 对, 前序未匹配: {n_prev - len(used_prev)}, 当前未匹配: {n_curr - len(used_curr)}")
    return matches


def track_eddies(eddy_series, threshold=base_threshold, max_daily_distance=max_daily_distance,
                 split_threshold=split_threshold, merge_threshold=merge_threshold):
    """
    全局涡旋轨迹追踪（基于多时间步匹配结果，支持轨迹中断续接）
    参数：
        eddy_series (list[list[dict]]): 时间序列涡旋数据（load_eddy_series的输出，已转换经度）
        threshold (float): 直接匹配阈值（默认全局base_threshold）
        max_daily_distance (float): 日均最大移动距离（默认全局max_daily_distance）
        split_threshold (float): 分裂阈值（默认全局split_threshold）
        merge_threshold (float): 合并阈值（默认全局merge_threshold）
    返回：
        defaultdict[list] | None: 轨迹字典（key=轨迹ID，value=轨迹点列表）；数据不足返回None
    轨迹追踪逻辑：
        1. 初始化第0时间步涡旋为初始轨迹；
        2. 逐时间步调用match_eddies_with_split_merge进行匹配，更新轨迹；
        3. 未匹配轨迹进入中断宽容期，若在宽容期内重新匹配则续接，否则终止；
        4. 未匹配当前涡旋作为新轨迹起始；
    轨迹点规范：
        每个轨迹点含11个必选字段（确保后续统计/可视化兼容性）：
        time_step/obs_id/latitude/longitude/amplitude/effective_radius/effective_area/date/event/contour_lat/contour_lon
    """
    # 轨迹点必选字段（确保所有轨迹点格式统一，避免后续处理报错）
    REQUIRED_FIELDS = {
        'time_step', 'obs_id', 'latitude', 'longitude', 'amplitude',
        'effective_radius', 'effective_area', 'date', 'event',
        'contour_lat', 'contour_lon'
    }

    # 至少需要2个时间步才能追踪轨迹（1个时间步无连续关系）
    if len(eddy_series) < 2:
        logger.error("轨迹追踪需要至少2个时间步的涡旋数据，当前数据不足")
        return None

    logger.info(f"开始全局涡旋轨迹追踪（日均最大移动距离: {max_daily_distance}km，直接匹配阈值: {threshold}）")

    # 初始化轨迹存储结构
    trajectories = defaultdict(list)  # 轨迹字典（key=轨迹ID，value=轨迹点列表）
    first_eddies = eddy_series[0]     # 第0时间步涡旋（初始轨迹来源）
    traj_counter = 0                  # 轨迹ID计数器（用于生成唯一ID）
    prev_to_traj = {}                 # 前序时间步涡旋索引→轨迹ID的映射（用于快速匹配）
    inactive_trajectories = defaultdict(int)  # 中断轨迹计数器（key=轨迹ID，value=中断步数）

    # --------------------------
    # 步骤1：初始化第0时间步轨迹
    # --------------------------
    for eddy_idx, eddy in enumerate(first_eddies):
        # 生成唯一轨迹ID（格式：traj_0, traj_1, ...）
        traj_id = f"traj_{traj_counter}"
        traj_counter += 1

        # 构建初始轨迹点（按必选字段填充）
        init_point = {
            'time_step': eddy['time_step'],  # 全局时间步（0）
            'obs_id': eddy['obs_id'],        # 当前文件内涡旋索引
            'latitude': eddy['contour_lat_mean'],  # 轮廓纬度均值（优先于中心点）
            'longitude': eddy['contour_lon_mean'], # 轮廓经度均值
            'amplitude': eddy['amplitude'],  # 振幅
            'effective_radius': eddy['effective_radius'],  # 有效半径
            'effective_area': eddy['effective_area'],      # 有效面积
            'date': eddy['date'],            # 日期字符串
            'event': 'start',                # 事件类型（初始轨迹标记为"start"）
            'contour_lat': eddy['contour_latitude'],  # 完整轮廓纬度（None若无）
            'contour_lon': eddy['contour_longitude']   # 完整轮廓经度
        }

        # 检查并补充缺失字段（避免NC文件变量缺失导致的KeyError）
        missing_fields = REQUIRED_FIELDS - set(init_point.keys())
        if missing_fields:
            logger.warning(f"初始轨迹点（轨迹ID: {traj_id}）缺失字段: {missing_fields}，自动补充默认值")
            for field in missing_fields:
                if field == 'effective_area':
                    # 面积缺失时从涡旋数据补充（若仍缺失则设为0.0）
                    init_point[field] = eddy.get('effective_area', 0.0)
                else:
                    # 其他字段缺失设为None（后续处理需兼容）
                    init_point[field] = None

        # 添加初始轨迹点到轨迹字典
        trajectories[traj_id].append(init_point)
        # 建立前序涡旋索引→轨迹ID的映射（第0时间步无前置，映射自身索引）
        prev_to_traj[eddy_idx] = traj_id

    logger.info(f"初始化完成：第0时间步含 {len(first_eddies)} 个涡旋，创建 {len(trajectories)} 条初始轨迹")

    # --------------------------
    # 步骤2：逐时间步追踪轨迹（从第1时间步开始）
    # --------------------------
    for time_step in tqdm(range(1, len(eddy_series)), desc="时间步轨迹追踪进度"):
        prev_eddies = eddy_series[time_step - 1]  # 前序时间步涡旋（t-1）
        curr_eddies = eddy_series[time_step]      # 当前时间步涡旋（t）
        curr_to_traj = {}                         # 当前时间步涡旋索引→轨迹ID的映射

        # 提取前序与当前时间步日期（用于日志）
        prev_date = prev_eddies[0]['date'] if prev_eddies else "unknown"
        curr_date = curr_eddies[0]['date'] if curr_eddies else "unknown"
        logger.info(f"\n开始匹配时间步 {time_step - 1} ({prev_date}) → {time_step} ({curr_date})")
        logger.info(f"前序时间步涡旋数: {len(prev_eddies)}, 当前时间步涡旋数: {len(curr_eddies)}")

        # 调用单时间步匹配函数，获取匹配结果
        matches = match_eddies_with_split_merge(
            prev_eddies, curr_eddies,
            threshold=threshold,
            max_daily_distance=max_daily_distance,
            split_threshold=split_threshold,
            merge_threshold=merge_threshold
        )

        # --------------------------
        # 步骤2.1：处理匹配结果（更新轨迹）
        # --------------------------
        for prev_idx, curr_idx, match_type in matches:
            # 前序涡旋对应的轨迹ID（从prev_to_traj映射获取）
            if prev_idx not in prev_to_traj:
                logger.warning(f"前序涡旋索引 {prev_idx} 无对应轨迹ID，跳过该匹配对")
                continue

            original_traj_id = prev_to_traj[prev_idx]
            curr_eddy = curr_eddies[curr_idx]  # 当前匹配的涡旋数据

            # 情况1：分裂事件（前序1个轨迹→当前多个轨迹，需创建新轨迹）
            if match_type == 'split':
                # 生成新轨迹ID（分裂分支轨迹）
                new_traj_id = f"traj_{traj_counter}"
                traj_counter += 1

                # 复制原轨迹的历史点（分裂轨迹继承原轨迹的所有前置点）
                trajectories[new_traj_id] = [p.copy() for p in trajectories[original_traj_id]]

                # 构建分裂事件的轨迹点
                split_point = {
                    'time_step': curr_eddy['time_step'],
                    'obs_id': curr_eddy['obs_id'],
                    'latitude': curr_eddy['contour_lat_mean'],
                    'longitude': curr_eddy['contour_lon_mean'],
                    'amplitude': curr_eddy['amplitude'],
                    'effective_radius': curr_eddy['effective_radius'],
                    'effective_area': curr_eddy['effective_area'],
                    'date': curr_eddy['date'],
                    'event': 'split',  # 事件类型标记为"split"
                    'contour_lat': curr_eddy['contour_latitude'],
                    'contour_lon': curr_eddy['contour_longitude']
                }

                # 补充缺失字段
                missing_fields = REQUIRED_FIELDS - set(split_point.keys())
                if missing_fields:
                    logger.warning(f"分裂轨迹点（轨迹ID: {new_traj_id}）缺失字段: {missing_fields}，自动补充")
                    for field in missing_fields:
                        if field == 'effective_area':
                            split_point[field] = curr_eddy.get('effective_area', 0.0)
                        else:
                            split_point[field] = None

                # 添加分裂点到新轨迹
                trajectories[new_traj_id].append(split_point)
                # 建立当前涡旋索引→新轨迹ID的映射
                curr_to_traj[curr_idx] = new_traj_id
                # 分裂轨迹为活跃状态，从中断列表移除（若存在）
                inactive_trajectories.pop(new_traj_id, None)

            # 情况2：合并事件（前序多个轨迹→当前1个轨迹，合并到其中一个轨迹）
            elif match_type == 'merge':
                # 若当前涡旋未被匹配，将其关联到原轨迹（避免重复合并）
                if curr_idx not in curr_to_traj:
                    curr_to_traj[curr_idx] = original_traj_id

                    # 构建合并事件的轨迹点
                    merge_point = {
                        'time_step': curr_eddy['time_step'],
                        'obs_id': curr_eddy['obs_id'],
                        'latitude': curr_eddy['contour_lat_mean'],
                        'longitude': curr_eddy['contour_lon_mean'],
                        'amplitude': curr_eddy['amplitude'],
                        'effective_radius': curr_eddy['effective_radius'],
                        'effective_area': curr_eddy['effective_area'],
                        'date': curr_eddy['date'],
                        'event': 'merge',  # 事件类型标记为"merge"
                        'contour_lat': curr_eddy['contour_latitude'],
                        'contour_lon': curr_eddy['contour_longitude']
                    }

                    # 补充缺失字段
                    missing_fields = REQUIRED_FIELDS - set(merge_point.keys())
                    if missing_fields:
                        logger.warning(f"合并轨迹点（轨迹ID: {original_traj_id}）缺失字段: {missing_fields}，自动补充")
                        for field in missing_fields:
                            if field == 'effective_area':
                                merge_point[field] = curr_eddy.get('effective_area', 0.0)
                            else:
                                merge_point[field] = None

                    # 添加合并点到原轨迹
                    trajectories[original_traj_id].append(merge_point)
                else:
                    # 若当前涡旋已被匹配（多前序合并到同一当前涡旋），标记为"merged"
                    trajectories[original_traj_id][-1]['event'] = 'merged'

                # 合并轨迹为活跃状态，从中断列表移除
                inactive_trajectories.pop(original_traj_id, None)

            # 情况3：直接匹配（前序1个→当前1个，延续原轨迹）
            else:
                # 构建直接匹配的轨迹点
                direct_point = {
                    'time_step': curr_eddy['time_step'],
                    'obs_id': curr_eddy['obs_id'],
                    'latitude': curr_eddy['contour_lat_mean'],
                    'longitude': curr_eddy['contour_lon_mean'],
                    'amplitude': curr_eddy['amplitude'],
                    'effective_radius': curr_eddy['effective_radius'],
                    'effective_area': curr_eddy['effective_area'],
                    'date': curr_eddy['date'],
                    'event': 'continue',  # 事件类型标记为"continue"
                    'contour_lat': curr_eddy['contour_latitude'],
                    'contour_lon': curr_eddy['contour_longitude']
                }

                # 补充缺失字段
                missing_fields = REQUIRED_FIELDS - set(direct_point.keys())
                if missing_fields:
                    logger.warning(f"直接匹配轨迹点（轨迹ID: {original_traj_id}）缺失字段: {missing_fields}，自动补充")
                    for field in missing_fields:
                        if field == 'effective_area':
                            direct_point[field] = curr_eddy.get('effective_area', 0.0)
                        else:
                            direct_point[field] = None

                # 添加直接匹配点到原轨迹
                trajectories[original_traj_id].append(direct_point)
                # 建立当前涡旋索引→原轨迹ID的映射
                curr_to_traj[curr_idx] = original_traj_id
                # 直接匹配轨迹为活跃状态，从中断列表移除
                inactive_trajectories.pop(original_traj_id, None)

        # --------------------------
        # 步骤2.2：处理未匹配轨迹（进入中断宽容期）
        # --------------------------
        for prev_idx, traj_id in prev_to_traj.items():
            # 前序涡旋未被匹配（无对应的当前涡旋），轨迹进入中断状态
            if prev_idx not in [m[0] for m in matches]:
                inactive_trajectories[traj_id] += 1  # 中断步数+1
                logger.debug(f"轨迹 {traj_id} 中断，当前中断步数: {inactive_trajectories[traj_id]}")

        # --------------------------
        # 步骤2.3：中断轨迹续接（宽容期内重新匹配）
        # --------------------------
        if inactive_trajectories and curr_eddies:
            # 提取宽容期内的中断轨迹特征（中断步数≤tolerance_steps）
            inactive_features = []  # 中断轨迹最后一个点的特征
            inactive_ids = []       # 中断轨迹ID列表
            for traj_id, inactive_steps in inactive_trajectories.items():
                if inactive_steps <= tolerance_steps:
                    # 获取轨迹最后一个点（中断前的最新点）
                    last_point = trajectories[traj_id][-1]
                    # 补充有效面积（避免后续相似度计算报错）
                    if 'effective_area' not in last_point:
                        logger.warning(f"中断轨迹 {traj_id} 的最后一个点缺失effective_area，补充为0.0")
                        last_point['effective_area'] = 0.0

                    # 提取特征用于相似度计算
                    inactive_features.append({
                        'lat': last_point['latitude'],
                        'lon': last_point['longitude'],
                        'radius': last_point['effective_radius'],
                        'amp': last_point['amplitude'],
                        'area': last_point['effective_area']
                    })
                    inactive_ids.append(traj_id)

            # 提取当前时间步未匹配涡旋的特征
            unmatched_curr_indices = [j for j in range(len(curr_eddies)) if j not in curr_to_traj]
            if inactive_features and unmatched_curr_indices:
                curr_features = [
                    {
                        'lat': curr_eddies[j]['contour_lat_mean'],
                        'lon': curr_eddies[j]['contour_lon_mean'],
                        'radius': curr_eddies[j]['effective_radius'],
                        'amp': curr_eddies[j]['amplitude'],
                        'area': curr_eddies[j]['effective_area'],
                        'idx': j  # 当前涡旋索引
                    }
                    for j in unmatched_curr_indices
                ]

                # 转换特征为NumPy数组（适配Numba加速函数）
                lats1 = np.array([f['lat'] for f in inactive_features])
                lons1 = np.array([f['lon'] for f in inactive_features])
                radii1 = np.array([f['radius'] for f in inactive_features])
                amps1 = np.array([f['amp'] for f in inactive_features])
                areas1 = np.array([f['area'] for f in inactive_features])

                lats2 = np.array([f['lat'] for f in curr_features])
                lons2 = np.array([f['lon'] for f in curr_features])
                radii2 = np.array([f['radius'] for f in curr_features])
                amps2 = np.array([f['amp'] for f in curr_features])
                areas2 = np.array([f['area'] for f in curr_features])

                # 计算中断轨迹与未匹配当前涡旋的相似度（放宽约束：距离×(宽容步数+1)，阈值×0.8）
                tolerance_similarity = calculate_similarity_matrix(
                    lats1, lons1, radii1, amps1, areas1,
                    lats2, lons2, radii2, amps2, areas2,
                    max_daily_distance=max_daily_distance * (tolerance_steps + 1)
                )

                # 按相似度降序匹配，续接轨迹
                flat_tol_indices = np.argsort(tolerance_similarity.flatten())[::-1]
                used_inactive = set()  # 已匹配的中断轨迹索引
                used_curr_tol = set()  # 已匹配的当前涡旋索引

                for idx in flat_tol_indices:
                    i = idx // len(curr_features)  # 中断轨迹索引
                    j = idx % len(curr_features)   # 当前涡旋特征索引

                    # 相似度超过放宽阈值（原阈值×0.8），且均未使用
                    if i not in used_inactive and j not in used_curr_tol and tolerance_similarity[i, j] > (threshold * 0.8):
                        traj_id = inactive_ids[i]          # 中断轨迹ID
                        curr_j = curr_features[j]['idx']   # 当前涡旋索引
                        curr_eddy = curr_eddies[curr_j]    # 当前涡旋数据

                        # 构建续接轨迹点（事件类型标记为"resumed"）
                        resumed_point = {
                            'time_step': curr_eddy['time_step'],
                            'obs_id': curr_eddy['obs_id'],
                            'latitude': curr_eddy['contour_lat_mean'],
                            'longitude': curr_eddy['contour_lon_mean'],
                            'amplitude': curr_eddy['amplitude'],
                            'effective_radius': curr_eddy['effective_radius'],
                            'effective_area': curr_eddy['effective_area'],
                            'date': curr_eddy['date'],
                            'event': 'resumed',  # 续接事件标记
                            'contour_lat': curr_eddy['contour_latitude'],
                            'contour_lon': curr_eddy['contour_longitude']
                        }

                        # 补充缺失字段
                        missing_fields = REQUIRED_FIELDS - set(resumed_point.keys())
                        if missing_fields:
                            logger.warning(f"续接轨迹点（轨迹ID: {traj_id}）缺失字段: {missing_fields}，自动补充")
                            for field in missing_fields:
                                if field == 'effective_area':
                                    resumed_point[field] = curr_eddy.get('effective_area', 0.0)
                                else:
                                    resumed_point[field] = None

                        # 添加续接点点到轨迹
                        trajectories[traj_id].append(resumed_point)
                        # 建立当前涡旋索引→轨迹ID的映射
                        curr_to_traj[curr_j] = traj_id
                        # 标记为已使用，避免重复匹配
                        used_inactive.add(i)
                        used_curr_tol.add(j)
                        # 续接后轨迹活跃，从中断列表移除
                        inactive_trajectories.pop(traj_id)

            logger.info(f"时间步 {time_step} 续接中断轨迹: {len(used_inactive)} 条")

        # --------------------------
        # 步骤2.4：处理未匹配当前涡旋（创建新轨迹）
        # --------------------------
        new_traj_count = 0  # 新增轨迹计数器
        for curr_idx in range(len(curr_eddies)):
            if curr_idx not in curr_to_traj:
                # 当前涡旋未被匹配（既不是直接匹配/分裂/合并，也不是中断续接），创建新轨迹
                curr_eddy = curr_eddies[curr_idx]
                new_traj_id = f"traj_{traj_counter}"
                traj_counter += 1

                # 构建新轨迹的初始点（事件类型标记为"start"）
                new_point = {
                    'time_step': curr_eddy['time_step'],
                    'obs_id': curr_eddy['obs_id'],
                    'latitude': curr_eddy['contour_lat_mean'],
                    'longitude': curr_eddy['contour_lon_mean'],
                    'amplitude': curr_eddy['amplitude'],
                    'effective_radius': curr_eddy['effective_radius'],
                    'effective_area': curr_eddy['effective_area'],
                    'date': curr_eddy['date'],
                    'event': 'start',  # 新轨迹标记
                    'contour_lat': curr_eddy['contour_latitude'],
                    'contour_lon': curr_eddy['contour_longitude']
                }

                # 补充缺失字段
                missing_fields = REQUIRED_FIELDS - set(new_point.keys())
                if missing_fields:
                    logger.warning(f"新增轨迹点（轨迹ID: {new_traj_id}）缺失字段: {missing_fields}，自动补充")
                    for field in missing_fields:
                        if field == 'effective_area':
                            new_point[field] = curr_eddy.get('effective_area', 0.0)
                        else:
                            new_point[field] = None

                # 添加新轨迹点到轨迹字典
                trajectories[new_traj_id].append(new_point)
                # 建立当前涡旋索引→新轨迹ID的映射
                curr_to_traj[curr_idx] = new_traj_id
                new_traj_count += 1

        # --------------------------
        # 步骤2.5：清理超宽容期的中断轨迹（终止）
        # --------------------------
        # 中断步数超过tolerance_steps的轨迹，视为永久终止，从轨迹字典移除
        to_remove = [traj_id for traj_id, steps in inactive_trajectories.items() if steps > tolerance_steps]
        for traj_id in to_remove:
            del inactive_trajectories[traj_id]
            # 移除轨迹前记录日志（可选：保留短轨迹需注释此句）
            # del trajectories[traj_id]  # 若需保留所有轨迹，注释此句；仅保留有效轨迹则保留
            logger.debug(f"轨迹 {traj_id} 中断超过宽容期（{tolerance_steps}步），已终止")

        # 更新前序涡旋→轨迹ID的映射（为下一时间步准备）
        prev_to_traj = curr_to_traj
        # 记录当前时间步轨迹统计
        logger.info(
            f"时间步 {time_step} 轨迹统计: 匹配 {len(matches)} 条, 续接 {len(used_inactive)} 条, 新增 {new_traj_count} 条, 总轨迹数: {len(trajectories)}")

    # --------------------------
    # 步骤3：最终轨迹点完整性检查（避免后续处理报错）
    # --------------------------
    for traj_id, traj_points in trajectories.items():
        for point_idx, point in enumerate(traj_points):
            # 强制补充effective_area（关键字段，缺失会导致统计/可视化错误）
            if 'effective_area' not in point:
                logger.error(f"轨迹 {traj_id} 的第 {point_idx} 个点缺失effective_area，强制补充为0.0")
                point['effective_area'] = 0.0
            # 检查必选字段数量（确保11个字段完整）
            if len(point) != len(REQUIRED_FIELDS):
                logger.warning(f"轨迹 {traj_id} 的第 {point_idx} 个点字段数量异常（{len(point)}/{len(REQUIRED_FIELDS)}），但已包含所有必选字段")

    logger.info(f"全局轨迹追踪完成，共生成 {len(trajectories)} 条轨迹")
    return trajectories

# --------------------------
# 8. 轨迹统计分析模块（输出关键统计指标）
# --------------------------
def analyze_trajectory_lengths(trajectories):
    """
    分析轨迹长度分布（时间步数），统计长轨迹占比、分裂/合并事件数量
    参数：
        trajectories (defaultdict[list]): 轨迹字典（track_eddies的输出）
    返回：
        Counter | None: 轨迹长度计数器（key=轨迹长度，value=对应长度的轨迹数）；无轨迹返回None
    统计指标：
        1. 总轨迹数、平均长度、最长轨迹长度；
        2. 长轨迹阈值（前25%分位数）及占比；
        3. 分裂/合并事件总数量；
    """
    if not trajectories:
        logger.warning("无轨迹数据可分析，返回空结果")
        return None

    # 提取所有轨迹的长度（时间步数）
    traj_lengths = [len(traj) for traj in trajectories.values()]
    # 统计各长度的轨迹数量（Counter自动计数）
    length_counts = Counter(traj_lengths)
    total_traj = len(trajectories)  # 总轨迹数
    # 定义长轨迹阈值（前25%分位数，即长度≥该值的轨迹为长轨迹）
    long_traj_threshold = np.percentile(traj_lengths, 75)

    # 输出统计日志（按规范格式）
    logger.info("\n===== 涡旋轨迹长度统计报告 =====")
    logger.info(f"总轨迹数: {total_traj}")
    logger.info(f"长轨迹阈值（前25%分位数）: {long_traj_threshold:.1f} 时间步")
    logger.info("轨迹长度 → 轨迹数量 (占比)")
    # 按长度升序输出统计
    for length in sorted(length_counts.keys()):
        count = length_counts[length]
        percentage = 100 * count / total_traj
        # 标记长轨迹（长度≥阈值）
        long_mark = '  ← 长轨迹' if length >= long_traj_threshold else ''
        logger.info(f"{length:>2d} 时间步 → {count:>4d} 条 ({percentage:.1f}%){long_mark}")

    # 计算关键统计量
    avg_length = np.mean(traj_lengths)  # 平均长度
    max_length = max(traj_lengths)      # 最长长度
    # 长轨迹数量与占比
    long_traj_count = sum(count for length, count in length_counts.items() if length >= long_traj_threshold)
    long_traj_ratio = 100 * long_traj_count / total_traj

    logger.info(f"\n===== 关键统计指标 =====")
    logger.info(f"平均轨迹长度: {avg_length:.2f} 时间步")
    logger.info(f"最长轨迹长度: {max_length} 时间步")
    logger.info(f"长轨迹数量: {long_traj_count} 条，占比: {long_traj_ratio:.1f}%")

    # 统计分裂/合并事件数量
    split_count = 0  # 分裂事件数
    merge_count = 0  # 合并事件数
    for traj in trajectories.values():
        for point in traj:
            event = point.get('event', '')
            if event == 'split':
                split_count += 1
            elif event == 'merge' or event == 'merged':
                merge_count += 1

    logger.info(f"分裂事件总数: {split_count} 次")
    logger.info(f"合并事件总数: {merge_count} 次")

    return length_counts

# --------------------------
# 9. 可视化模块（绘制轨迹图与时间步快照）
# --------------------------
def add_basemap(ax, basemap_path=None, lon_range=(-45, 0), lat_range=(0, 30)):
    """
    为绘图轴（ax）添加底图（支持NC格式ADT数据或默认地理底图）
    参数：
        ax (matplotlib.axes.Axes): 绘图轴对象
        basemap_path (str | None): NC格式底图文件目录（如ADT数据）；None则使用默认地理底图
        lon_range (tuple): 经度范围（min_lon, max_lon）
        lat_range (tuple): 纬度范围（min_lat, max_lat）
    返回：
        Basemap: 底图对象（用于后续绘图）
    底图逻辑：
        1. 若basemap_path为None：绘制默认地理底图（海岸线、经纬线、网格）；
        2. 若为目录：按调用次数读取目录内NC文件（ADT数据），绘制彩色等高线（排除陆地值）；
    关键优化：
        排除陆地值（-2147483647）和NaN，使用95%分位数确定色阶范围（避免异常值影响）
    """
    # 用函数属性记录调用次数（静态变量，实现NC文件按时间步顺序读取）
    if not hasattr(add_basemap, 'call_count'):
        add_basemap.call_count = 0  # 初始化调用计数器（从0开始）

    # 情况1：使用默认地理底图（无NC底图文件）
    if basemap_path is None:
        m = Basemap(
            projection='cyl',  # 圆柱投影（适合小范围区域，无变形）
            llcrnrlat=lat_range[0], urcrnrlat=lat_range[1],  # 纬度范围
            llcrnrlon=lon_range[0], urcrnrlon=lon_range[1],  # 经度范围
            resolution='l',  # 分辨率（l=低，m=中，h=高；低分辨率速度快）
            ax=ax
        )
        # 绘制基础地理元素
        m.drawcoastlines(linewidth=0.8, color='black')  # 海岸线
        m.drawparallels(np.arange(-90, 91, 15), labels=[1, 0, 0, 0], fontsize=8)  # 纬线（下侧标签）
        m.drawmeridians(np.arange(-180, 181, 30), labels=[0, 0, 0, 1], fontsize=8)  # 经线（右侧标签）
        ax.grid(True, linestyle='--', alpha=0.7)  # 网格线（虚线，透明度0.7）
        return m

    # 情况2：使用NC格式底图文件（如ADT数据）
    # 检查目录有效性
    if not os.path.isdir(basemap_path):
        logger.error(f"底图目录 {basemap_path} 不是有效路径，使用默认底图")
        return add_basemap(ax, None, lon_range, lat_range)

    # 筛选目录下所有NC文件并按名称排序（假设文件名含时间信息，如"adt_20230101.nc"）
    nc_files = sorted([f for f in os.listdir(basemap_path) if f.lower().endswith('.nc')])
    if not nc_files:
        logger.warning(f"底图目录 {basemap_path} 下未找到NC文件，使用默认底图")
        return add_basemap(ax, None, lon_range, lat_range)

    # 按调用次数选择NC文件（实现时间步与底图文件对应）
    current_call = add_basemap.call_count
    # 若调用次数超过文件数，使用最后一个文件
    if current_call >= len(nc_files):
        logger.error(f"底图调用次数({current_call + 1})超过NC文件数({len(nc_files)})，使用最后一个文件")
        current_call = len(nc_files) - 1

    # 读取当前NC文件
    nc_file_path = os.path.join(basemap_path, nc_files[current_call])
    try:
        with h5netcdf.File(nc_file_path, 'r') as ds:
            # 创建基础Basemap对象
            m = Basemap(
                projection='cyl',
                llcrnrlat=lat_range[0], urcrnrlat=lat_range[1],
                llcrnrlon=lon_range[0], urcrnrlon=lon_range[1],
                resolution='l',
                ax=ax
            )
            # 绘制基础地理元素
            m.drawcoastlines(linewidth=0.8, color='black')
            m.drawparallels(np.arange(-90, 91, 15), labels=[1, 0, 0, 0], fontsize=8)
            m.drawmeridians(np.arange(-180, 181, 30), labels=[0, 0, 0, 1], fontsize=8)

            # 提取ADT数据（尝试常见变量名，确保兼容性）
            lat = ds['latitude'][:]
            lon = ds['longitude'][:]
            adt_var_names = ['adt', 'sla', 'zos', 'ssh']  # 常见ADT/SLA变量名
            adt_data = None
            for var_name in adt_var_names:
                if var_name in ds.variables:
                    adt_data = ds[var_name][:].squeeze()  # 挤压维度（去除长度为1的维度）
                    break

            # 若未找到ADT变量，返回默认底图
            if adt_data is None:
                logger.warning(f"NC底图文件 {nc_files[current_call]} 中未找到ADT相关变量（尝试: {adt_var_names}），使用默认底图")
                return m

            # --------------------------
            # 关键优化：处理陆地值，避免影响色阶
            # --------------------------
            # 1. 创建陆地值掩码（NC文件中常用-2147483647表示陆地，同时排除NaN）
            land_mask = (adt_data == -2147483647) | np.isnan(adt_data)
            adt_masked = np.ma.masked_where(land_mask, adt_data)  # 掩码陆地值

            # 2. 计算有效海洋数据的统计量（仅用海洋数据确定色阶，避免陆地值干扰）
            valid_adt = adt_masked.compressed()  # 提取所有非掩码数据（仅海洋）
            if len(valid_adt) > 0:
                # 使用95%分位数确定色阶范围（排除5%异常值，使颜色分布更均匀）
                q05, q95 = np.percentile(valid_adt, [5, 95])
                adt_min, adt_max = q05, q95

                # 确保色阶范围有效（避免所有值相等导致色阶失效）
                if adt_min == adt_max:
                    adt_min -= 0.1
                    adt_max += 0.1

                # 扩展10%边距，提升颜色过渡效果
                margin = (adt_max - adt_min) * 0.1
                adt_min -= margin
                adt_max += margin

                # 生成色阶（40级，确保颜色过渡平滑）
                levels = np.linspace(adt_min, adt_max, 40)
            else:
                # 无有效海洋数据时，使用默认色阶
                levels = np.linspace(-1, 1, 40)
                logger.warning(f"NC底图文件 {nc_files[current_call]} 中无有效海洋ADT数据，使用默认色阶")

            # 3. 绘制ADT彩色等高线
            lon_grid, lat_grid = np.meshgrid(lon, lat)  # 生成经纬度网格
            im = m.contourf(
                lon_grid, lat_grid, adt_masked,
                levels=levels,
                cmap='coolwarm',  # 冷暖色映射（适合表示正负异常）
                alpha=0.8,        # 透明度（避免遮挡轨迹）
                extend='both',    # 超出色阶范围的数值用极值颜色表示
                antialiased=True  # 抗锯齿，使边缘更平滑
            )

            # 添加颜色条（底部放置，避免遮挡轨迹）
            cbar = m.colorbar(im, location='bottom', pad="10%")
            cbar.set_label('绝对动力地形 (m)', fontsize=8)  # 颜色条标签
            cbar.ax.tick_params(labelsize=6)  # 颜色条刻度字体大小

            # 记录底图加载信息
            logger.info(
                f"底图加载完成：第{current_call + 1}次调用，文件: {nc_files[current_call]}，ADT范围: [{adt_min:.4f}, {adt_max:.4f}] m")
            add_basemap.call_count += 1  # 更新调用次数（下一次调用读取下一个文件）
            return m

    except Exception as e:
        # 捕获NC文件读取异常，记录错误并返回默认底图
        logger.error(f"加载NC底图文件 {nc_files[current_call]} 失败: {str(e)}，使用默认底图", exc_info=True)
        add_basemap.call_count += 1  # 即使失败也更新调用次数，避免重复尝试
        return m


def _plot_nc_data(ax, ds):
    """
    内部辅助函数：从NC数据集提取并绘制地理数据（用于兼容多种NC格式）
    参数：
        ax (matplotlib.axes.Axes): 绘图轴对象
        ds (h5netcdf.Dataset): 已打开的NC数据集
    说明：
        1. 自动识别经纬度变量（含"lon"/"lat"关键字）；
        2. 选择第一个非经纬度的2D变量作为数据变量；
        3. 绘制数据为伪彩色图并添加颜色条
    """
    # 自动识别经纬度变量
    lon_var = None
    lat_var = None
    for var in ds.variables:
        if 'lon' in var.lower() and not lon_var:
            lon_var = var
        if 'lat' in var.lower() and not lat_var:
            lat_var = var

    # 自动识别数据变量（优先选择2D变量）
    data_var = None
    for var in ds.variables:
        if var != lon_var and var != lat_var and len(ds[var].shape) == 2:
            data_var = var
            break

    # 变量识别失败时返回
    if not all([lon_var, lat_var, data_var]):
        logger.warning("NC文件变量识别失败：未找到有效的经纬度或2D数据变量")
        return

    # 读取数据并处理维度顺序
    lons = ds[lon_var][:]
    lats = ds[lat_var][:]
    data = ds[data_var][:]
    # 确保数据维度为(lat, lon)，若不匹配则转置
    if data.shape != (len(lats), len(lons)):
        data = data.T

    # 绘制数据为伪彩色图
    im = ax.imshow(
        data,
        extent=[lons.min(), lons.max(), lats.min(), lats.max()],  # 坐标范围
        aspect='auto',  # 自动调整纵横比
        alpha=0.7,      # 透明度
        cmap='viridis'  # 通用颜色映射
    )

    # 添加颜色条
    cbar = plt.colorbar(im, ax=ax, orientation='vertical', fraction=0.046, pad=0.04)
    cbar.set_label(f'{data_var}')  # 颜色条标签为数据变量名


def draw_eddy_shape(ax, contour_lon, contour_lat, radius, color, alpha=0.3):
    """
    绘制涡旋形状（优先使用轮廓数据，无轮廓时用圆形近似）
    参数：
        ax (matplotlib.axes.Axes): 绘图轴对象
        contour_lon (np.ndarray | None): 涡旋轮廓经度（数组）
        contour_lat (np.ndarray | None): 涡旋轮廓纬度（数组）
        radius (float): 涡旋有效半径（km，用于圆形近似）
        color (tuple): 颜色RGB值（如(0.1, 0.2, 0.3, 1.0)）
        alpha (float): 透明度（0~1）
    说明：
        1. 若轮廓数据有效（长度>1），绘制多边形轮廓；
        2. 否则，用有效半径绘制圆形（1°≈111km，转换半径为度）
    """
    # 情况1：使用轮廓数据绘制多边形
    if contour_lon is not None and contour_lat is not None and len(contour_lon) > 1:
        try:
            ax.fill(
                contour_lon, contour_lat,
                color=color,
                alpha=alpha,
                edgecolor='black',  # 轮廓边缘黑色
                linewidth=0.5       # 边缘线宽
            )
            return
        except Exception as e:
            logger.debug(f"轮廓绘制失败，使用圆形替代: {str(e)}")  # 轮廓数据异常时降级为圆形

    # 情况2：使用有效半径绘制圆形（近似）
    if radius > 0:
        # 转换半径单位：km → 度（1度≈111km）
        radius_deg = radius / 111.0
        # 计算中心坐标（轮廓均值或0）
        center_lon = np.mean(contour_lon) if contour_lon is not None else 0
        center_lat = np.mean(contour_lat) if contour_lat is not None else 0
        # 创建圆形补丁并添加到轴
        circle = plt.Circle(
            (center_lon, center_lat),
            radius_deg,
            color=color,
            alpha=alpha,
            edgecolor='black',
            linewidth=0.5
        )
        ax.add_patch(circle)


def create_custom_legend(ax):
    """
    为轨迹图创建自定义图例（仅包含轨迹线和涡旋形状）
    参数：
        ax (matplotlib.axes.Axes): 绘图轴对象
    """
    # 定义图例元素（轨迹线和涡旋轮廓）
    legend_elements = [
        Line2D(
            [0], [0],
            color='gray',
            lw=1.5,
            alpha=0.8,
            label='轨迹线'
        ),
        Line2D(
            [0], [0],
            color='blue',
            lw=0,
            marker='s',
            markersize=8,
            markerfacecolor='blue',
            alpha=0.3,
            markeredgecolor='black',
            label='涡旋轮廓'
        )
    ]
    # 添加图例到轴
    ax.legend(
        handles=legend_elements,
        loc='upper right',
        fontsize=8,
        title="图例",
        title_fontsize=9,
        frameon=True,
        framealpha=0.9  # 图例背景半透明
    )


def visualize_trajectories(trajectories, eddy_series, output_dir='.', lon_range=None, lat_range=None, basemap_path=None,
                           max_pixels=60000):
    """
    绘制总轨迹图（所有轨迹叠加）并调用时间步可视化
    参数：
        trajectories (defaultdict[list]): 轨迹字典
        eddy_series (list[list[dict]]): 时间序列涡旋数据
        output_dir (str): 输出目录（默认当前目录）
        lon_range (tuple | None): 经度范围（min_lon, max_lon）；None则自动计算
        lat_range (tuple | None): 纬度范围（min_lat, max_lat）；None则自动计算
        basemap_path (str | None): 底图NC文件目录
        max_pixels (int): 最大像素限制（避免图像过大，默认60000）
    输出：
        保存总轨迹图到output_dir/eddy_trajectories_total.png
        调用visualize_time_steps生成时间步快照
    """
    # 筛选有效轨迹（长度>1，避免单点轨迹）
    valid_trajectories = {tid: traj for tid, traj in trajectories.items() if len(traj) > 1}
    if not valid_trajectories:
        logger.warning("无有效轨迹（长度>1）可可视化，跳过绘图")
        return

    # 自动计算坐标范围（若未指定）
    if lon_range is None or lat_range is None:
        all_lats = []
        all_lons = []
        for traj in valid_trajectories.values():
            for p in traj:
                all_lats.append(p['latitude'])
                all_lons.append(p['longitude'])
        # 使用1%和99%分位数排除极端值
        lon_percentile = np.percentile(all_lons, [1, 99])
        lat_percentile = np.percentile(all_lats, [1, 99])
        lon_range = [lon_percentile[0], lon_percentile[1]] if lon_range is None else lon_range
        lat_range = [lat_percentile[0], lat_percentile[1]] if lat_range is None else lat_range

    # 筛选可视化区域内的轨迹（至少有一个点在范围内）
    visual_traj = {}
    for tid, traj in valid_trajectories.items():
        in_range = False
        for p in traj:
            if (lon_range[0] <= p['longitude'] <= lon_range[1] and
                lat_range[0] <= p['latitude'] <= lat_range[1]):
                in_range = True
                break
        if in_range:
            visual_traj[tid] = traj

    logger.info(f"可视化区域内的有效轨迹数: {len(visual_traj)}")
    if not visual_traj:
        logger.warning("可视化区域内无有效轨迹，跳过绘图")
        return

    # 生成离散颜色映射（避免颜色过多导致混淆）
    num_colors = len(visual_traj)
    if num_colors <= 10:
        # 少量轨迹用鲜明的tab10颜色
        colors = plt.cm.tab10(np.linspace(0, 1, num_colors))
    else:
        # 大量轨迹用主色+亮度变化（循环使用10种主色，调整亮度区分）
        base_colors = plt.cm.tab10(np.linspace(0, 1, 10))
        colors = []
        for i in range(num_colors):
            base_idx = i % 10
            lightness_factor = (i // 10) % 5  # 0-4循环
            lightness = 0.8 + lightness_factor * 0.05  # 0.8-1.0亮度
            color = base_colors[base_idx] * lightness
            color = np.clip(color, 0.0, 1.0)  # 确保RGB值在0-1范围内
            color[-1] = 1.0  # 保持alpha通道为1
            colors.append(color)

    # 建立轨迹ID→颜色的映射
    traj_color_map = {tid: colors[i] for i, tid in enumerate(visual_traj.keys())}

    # 计算图像尺寸（限制最大像素）
    lon_span = lon_range[1] - lon_range[0]
    lat_span = lat_range[1] - lat_range[0]
    aspect_ratio = lon_span / lat_span if lat_span != 0 else 1.5
    fig_width = 12
    fig_height = fig_width / aspect_ratio
    # 限制高度，避免图像过大
    if fig_height * 100 > max_pixels:
        fig_height = max_pixels / 100
        fig_width = fig_height * aspect_ratio

    # 绘制总轨迹图
    fig = plt.figure(figsize=(fig_width, fig_height))
    ax = fig.add_subplot(111)
    # 添加底图
    add_basemap(ax, basemap_path, lon_range, lat_range)
    ax.set_xlim(lon_range)
    ax.set_ylim(lat_range)

    # 绘制每条轨迹的线和涡旋轮廓
    for tid, traj in visual_traj.items():
        color = traj_color_map[tid]
        # 轨迹点坐标
        lats = [p['latitude'] for p in traj]
        lons = [p['longitude'] for p in traj]
        # 绘制轨迹线
        ax.plot(lons, lats, '-', color=color, alpha=0.8, linewidth=1.5)
        # 绘制每个时间步的涡旋形状
        for p in traj:
            draw_eddy_shape(
                ax,
                p['contour_lon'],
                p['contour_lat'],
                p['effective_radius'],
                color,
                alpha=0.3
            )

    # 添加图例、标题和标签
    create_custom_legend(ax)
    ax.set_title('涡旋追踪总轨迹图（含涡旋轮廓）')
    ax.set_xlabel('经度')
    ax.set_ylabel('纬度')

    # 保存总轨迹图
    os.makedirs(output_dir, exist_ok=True)
    total_map_path = os.path.join(output_dir, 'eddy_trajectories_total.png')
    plt.savefig(total_map_path, dpi=200, bbox_inches='tight')
    logger.info(f"总轨迹图已保存至: {total_map_path}")

    # 调用时间步可视化（传递颜色映射和筛选后的轨迹）
    visualize_time_steps(
        trajectories, eddy_series, output_dir,
        lon_range, lat_range, basemap_path,
        traj_color_map, visual_traj
    )

    plt.close(fig)


def visualize_time_steps(trajectories, eddy_series, output_dir='.', lon_range=None, lat_range=None,
                         basemap_path=None, traj_color_map=None, visual_traj=None):
    """
    生成每个时间步的轨迹快照（展示轨迹随时间的演化）
    参数：
        trajectories (defaultdict[list]): 轨迹字典
        eddy_series (list[list[dict]]): 时间序列涡旋数据
        output_dir (str): 输出目录
        lon_range (tuple): 经度范围
        lat_range (tuple): 纬度范围
        basemap_path (str | None): 底图NC文件目录
        traj_color_map (dict): 轨迹ID→颜色的映射
        visual_traj (dict): 可视化区域内的轨迹
    输出：
        保存每个时间步的快照到output_dir/time_steps/time_step_*.png
    """
    os.makedirs(output_dir, exist_ok=True)
    time_steps = range(len(eddy_series))
    # 默认可视化所有有效轨迹（长度>1）
    if not visual_traj:
        visual_traj = {tid: traj for tid, traj in trajectories.items() if len(traj) > 1}

    # 创建时间步输出目录
    time_step_dir = os.path.join(output_dir, 'time_steps')
    os.makedirs(time_step_dir, exist_ok=True)

    # 确定坐标范围（若未指定）
    if lon_range is None or lat_range is None:
        all_lats = []
        all_lons = []
        for traj in visual_traj.values():
            for p in traj:
                all_lats.append(p['latitude'])
                all_lons.append(p['longitude'])
        lon_range = [min(all_lons), max(all_lons)] if lon_range is None else lon_range
        lat_range = [min(all_lats), max(all_lats)] if lat_range is None else lat_range

    # 计算图像尺寸
    lon_span = lon_range[1] - lon_range[0]
    lat_span = lat_range[1] - lat_range[0]
    aspect_ratio = lon_span / lat_span if lat_span != 0 else 1.5
    fig_width = 10
    fig_height = fig_width / aspect_ratio
    # 限制高度，避免过大
    if fig_height > 60:
        fig_height = 60
        fig_width = fig_height * aspect_ratio

    # 逐时间步生成快照
    for t in tqdm(time_steps, desc="生成时间步快照"):
        fig = plt.figure(figsize=(fig_width, fig_height))
        ax = fig.add_subplot(111)
        # 添加底图
        add_basemap(ax, basemap_path)
        ax.set_xlim(lon_range)
        ax.set_ylim(lat_range)

        # 绘制截至当前时间步的轨迹
        for tid, traj in visual_traj.items():
            color = traj_color_map.get(tid, (0.5, 0.5, 0.5, 1.0))  # 默认灰色
            # 筛选当前时间步及之前的轨迹点
            current_points = [p for p in traj if p['time_step'] <= t]
            if len(current_points) < 1:
                continue

            # 轨迹点坐标
            lats = [p['latitude'] for p in current_points]
            lons = [p['longitude'] for p in current_points]
            # 绘制轨迹线
            ax.plot(lons, lats, '-', color=color, alpha=0.8, linewidth=1.5)
            # 绘制涡旋轮廓
            for p in current_points:
                draw_eddy_shape(
                    ax,
                    p['contour_lon'],
                    p['contour_lat'],
                    p['effective_radius'],
                    color,
                    alpha=0.3
                )

        # 添加图例、标题和标签
        create_custom_legend(ax)
        # 获取当前时间步日期（若有）
        date_str = eddy_series[t][0]['date'] if (t < len(eddy_series) and eddy_series[t]) else "未知日期"
        ax.set_title(f'涡旋追踪结果 - 时间步 {t} ({date_str})')
        ax.set_xlabel('经度')
        ax.set_ylabel('纬度')

        # 保存当前时间步快照
        img_path = os.path.join(time_step_dir, f'time_step_{t}.png')
        plt.savefig(img_path, dpi=200, bbox_inches='tight')
        plt.close(fig)

    logger.info(f"所有时间步快照已保存至: {time_step_dir}")

# --------------------------
# 10. 结果保存模块（将轨迹数据写入NC文件）
# --------------------------
def save_trajectories_to_nc(trajectories, output_file):
    """
    将轨迹数据保存为NC文件（便于后续分析和共享）
    参数：
        trajectories (defaultdict[list]): 轨迹字典
        output_file (str): 输出NC文件路径（如"eddy_trajectories.nc"）
    说明：
        保存字段包括轨迹ID、长度、经纬度、振幅、半径、时间步、日期、事件类型
        缺失值用NaN或-1填充（兼容NC格式）
    """
    try:
        # 筛选有效轨迹（长度>1）
        valid_trajectories = {tid: traj for tid, traj in trajectories.items() if len(traj) > 1}
        n_trajectories = len(valid_trajectories)
        if n_trajectories == 0:
            logger.warning("无有效轨迹（长度>1）可保存，跳过写入")
            return

        # 确定最大轨迹长度（用于NC维度定义）
        max_length = max(len(traj) for traj in valid_trajectories.values())

        # 创建并写入NC文件
        with h5netcdf.File(output_file, 'w', decode_vlen_strings=False) as f:
            # 定义维度
            f.dimensions['trajectory'] = n_trajectories  # 轨迹数量
            f.dimensions['obs'] = max_length             # 最大观测数（时间步数）

            # 创建变量（按NC格式规范定义类型和填充值）
            f.create_variable('trajectory_id', ('trajectory',), 'S20')  # 轨迹ID（字符串）
            f.create_variable('trajectory_length', ('trajectory',), 'i4')  # 轨迹长度（整数）
            # 浮点变量用NaN填充，整数用-1填充
            f.create_variable('latitude', ('trajectory', 'obs'), 'f4', fillvalue=np.nan)
            f.create_variable('longitude', ('trajectory', 'obs'), 'f4', fillvalue=np.nan)
            f.create_variable('amplitude', ('trajectory', 'obs'), 'f4', fillvalue=np.nan)
            f.create_variable('effective_radius', ('trajectory', 'obs'), 'f4', fillvalue=np.nan)
            f.create_variable('time_step', ('trajectory', 'obs'), 'i4', fillvalue=-1)
            f.create_variable('date', ('trajectory', 'obs'), 'S10', fillvalue=b'unknown')  # 日期字符串
            f.create_variable('event', ('trajectory', 'obs'), 'S10', fillvalue=b'none')    # 事件类型

            # 填充数据
            for traj_idx, (tid, traj_points) in enumerate(valid_trajectories.items()):
                traj_len = len(traj_points)
                # 轨迹元信息
                f['trajectory_id'][traj_idx] = tid.encode('utf-8')  # 字符串需编码为字节
                f['trajectory_length'][traj_idx] = traj_len

                # 逐个轨迹点填充
                for point_idx, point in enumerate(traj_points):
                    f['latitude'][traj_idx, point_idx] = point['latitude']
                    f['longitude'][traj_idx, point_idx] = point['longitude']
                    f['amplitude'][traj_idx, point_idx] = point['amplitude']
                    f['effective_radius'][traj_idx, point_idx] = point['effective_radius']
                    f['time_step'][traj_idx, point_idx] = point['time_step']
                    f['date'][traj_idx, point_idx] = point['date'].encode('utf-8')
                    f['event'][traj_idx, point_idx] = point.get('event', 'none').encode('utf-8')

            # 添加全局属性（元数据）
            f.attrs['title'] = 'Eddy trajectories with split and merge events'
            f.attrs['history'] = f'Created on {datetime.now().strftime("%Y-%m-%d %H:%M:%S")}'
            f.attrs['max_daily_distance'] = f'{max_daily_distance}km (涡旋日均最大移动距离约束)'
            f.attrs['thresholds'] = f'base={base_threshold}, split={split_threshold}, merge={merge_threshold}'

        logger.info(f"轨迹数据已成功保存至NC文件: {output_file}，有效轨迹数: {n_trajectories}")

    except Exception as e:
        logger.error(f"保存轨迹至NC文件失败: {str(e)}", exc_info=True)

# --------------------------
# 11. 主程序（串联所有模块）
# --------------------------
def main():
    """
    主程序：涡旋轨迹追踪完整流程
    步骤：
    1. 配置参数（输入目录、输出路径、底图路径等）；
    2. 加载涡旋数据并转换经度范围；
    3. 执行轨迹追踪（含分裂/合并检测、中断续接）；
    4. 分析轨迹长度分布及事件统计；
    5. 保存轨迹数据至NC文件；
    6. 可视化轨迹（总轨迹图+时间步快照）
    """
    # 配置参数（根据实际数据路径修改）
    eddy_dir = "D:/EddyTrajData/EddyHex/2DC"  # 涡旋NC文件目录
    output_nc = "eddy_trajectories.nc"         # 输出轨迹NC文件
    output_dir = "eddy_tracking_results"       # 可视化结果输出目录
    basemap_path = "D:/EddyTrajData/EddyHex/SSH"  # 底图NC文件目录（ADT数据）

    # 日志输出配置信息
    logger.info("===== 涡旋轨迹追踪分析开始 =====")
    logger.info(f"涡旋数据目录: {eddy_dir}")
    logger.info(f"核心参数: 匹配阈值={base_threshold}, 日均最大移动距离={max_daily_distance}km")
    logger.info(f"事件阈值: 分裂={split_threshold}, 合并={merge_threshold}, 中断宽容期={tolerance_steps}步")

    try:
        # 步骤1：加载涡旋时间序列数据
        eddy_series = load_eddy_series(eddy_dir)
        if not eddy_series:
            logger.error("未加载到有效涡旋数据，程序终止")
            return

        # 步骤2：转换经度范围（360°→180°，适配可视化）
        eddy_series = check_and_convert_lon(eddy_series)

        # 步骤3：执行轨迹追踪
        trajectories = track_eddies(
            eddy_series,
            threshold=base_threshold,
            max_daily_distance=max_daily_distance,
            split_threshold=split_threshold,
            merge_threshold=merge_threshold
        )

        if not trajectories:
            logger.error("轨迹追踪失败，未生成任何轨迹")
            return

        # 步骤4：轨迹统计分析
        analyze_trajectory_lengths(trajectories)

        # 步骤5：保存轨迹数据
        save_trajectories_to_nc(trajectories, output_nc)

        # 步骤6：可视化轨迹（指定经纬度范围，可根据研究区域调整）
        visualize_trajectories(
            trajectories, eddy_series, output_dir,
            lon_range=[-45, 0], lat_range=[0, 30],  # 示例区域：大西洋部分区域
            basemap_path=basemap_path
        )

    except Exception as e:
        # 捕获全局异常，记录详细错误信息
        logger.error(f"程序执行过程中发生错误: {str(e)}", exc_info=True)

    logger.info("===== 涡旋轨迹追踪分析完成 =====")


# 程序入口（当脚本直接运行时执行主程序）
if __name__ == "__main__":
    main()